Spin-oscillator model for DNA/RNA unzipping by mechanical force 
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We model unzipping of DNA/RNA molecules subject to an external force by a spin-oscillator 
system. The system comprises a macroscopic degree of freedom, represented by a one-dimensional 
oscillator, and internal degrees of freedom, represented by Glauber spins with nearest-neighbor 
interaction and a coupling constant proportional to the oscillator position. At a critical value F c 
of an applied external force F, the oscillator rest position (order parameter) changes abruptly and 
the system undergoes a first-order phase transition. When the external force is cycled at different 
rates, the extension given by the oscillator position exhibits a hysteresis cycle at high loading rates 
whereas it moves reversibly over the equilibrium force-extension curve at very low loading rates. 
Under constant force, the logarithm of the residence time at the stable and metastable oscillator 
rest position is proportional to (F — F c ) as in an Arrhenius law. 

PACS numbers: 05.40.-a,05.50.+q,64.60.De,87.15.Cc 



I. INTRODUCTION 

Many physical situations can be modelled by a me- 
chanical system coupled to a thermal bath or to spin sys- 
tems. Examples abound, the collective Jahn- Teller effect 
has been analyzed by spin-phonon systems [IH3], mass 
spectrometry through a nanoelectromechanical oscillator 
whose resonant frequency decreases as single molecules 
are added thereto [4], decoherence of a spin representing 
a two-level system due to coupling to a boson bath (the 
spin-boson system) [5] , a classical oscillator coupled to a 
spin causes wave function collapse thereof (6], a 1/2-spin 
representing a nonlinear Josephson phase quantum bit is 
coupled to an oscillator (superconducting resonator) and 
to a classical signal [5] , rippling in clamped graphene 
sheets has been investigated by means of a spin-string 
system etc. 

Very recently, we have introduced a simple model in 
which a single oscillator is coupled to a chain of Ising 
spins undergoing Glauber dynamics in contact with a 
thermal bath [TUHT^] . In our model, the spins in the 
chain are coupled only to their nearest neighbors, but 
their coupling constant is proportional to the oscillator 
position, which makes their interaction effectively long 
ranged. In equilibrium, elimination of the oscillator co- 
ordinates gives rise to an effective spin interaction equiv- 
alent to a one dimensional Ising model with mean field 
coupling [13] . There is a second order phase transition 
at a finite temperature, with the oscillator rest position 
as its order parameter. Above the critical temperature, 
the oscillator rest position is zero, thereby coinciding 
with that of the uncoupled oscillator. Below the criti- 
cal temperature, two symmetric nonzero rest positions 
issue forth symmetrically from zero as in the case of a 
pitchfork bifurcation. In the limit of fast relaxation of 
the spins compared to the natural period of the oscilla- 



tor, the oscillator position satisfies an effective equation 
having both nonlinear force and nonlinear friction terms 
[TU1 IT2] . Interestingly, this nonlinear friction arises from 
the coupling of the macroscopic elastic mode with the 
internal degrees of freedom (modelled in our system by 
the spins). A related mechanism has been proposed to 
explain the "internal friction" observed in experiments 
with proteins or polymers in solution [14] . 

In recent years, technological development has allowed 
to manipulate or visualize individual molecules and to 
measure microscopic forces with high precision instru- 
ments. These single-molecule experiments (SME) pro- 
vide key information about the thermodynamic and ki- 
netic properties of biomolecules, offering a complemen- 
tary but different perspective to understand molecular 
processes. An extensive review of these techniques can 
be found in [T5]. Using SME, distributions describing 
certain molecular properties can be measured, thereby 
allowing to characterize the kinetics of biomolecular reac- 
tions and the observation of possible intermediate states. 
A typical outcome of SME are the force-extension curves 
for DNA, RNA and other biomolecules. In a seminal pa- 
per, Liphardt et al. [16] pull the nucleic acid molecule by 
an increasing applied force until it unfolds at a critical 
value of the force, F c ~ 14.5pN. Afterwards, the molecule 
is pushed back, by decreasing the force, until it refolds. 
At low pulling (pushing) rates, the stretching and re- 
laxing force-extension curves are superimposed, and the 
molecule unfolds at the critical value of the force F c . Tak- 
ing a closer look at this transition, hopping between the 
two possible extension values is observed. This suggests 
that the system is bistable: there are two possible states 
of the molecule with stochastic transitions between them. 
This physical picture is confirmed by experiments car- 
ried out at constant load, in a narrow region around the 
critical force (see fig. 2(C) of ref. [S]). When cycles 
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of pulling/pushing the molecule are carried out at high 
loading rates, the extension of the molecule occurs at a 
higher force F + > F c) whereas the hairpin folds at a lower 
force F_ < F c . Thus a hysteresis cycle arises, and some 
authors have been claimed this to be a signature of irre- 
versible non-equilibrium behavior |16H19j . More recently 
many works have tried to understand these unzipping 
experiments from a physical point of view |17H26| . 

In this paper, we add an external force to our previous 
oscillator-spin model |10) and analyze the resulting force- 
extension curves. Qualitatively, these curves have the 
same features as those of the force-extension curves mea- 
sured in experiments with DNA, RNA and other biopoly- 
mers described above. At subcritical temperatures, our 
spin-oscillator system has a first order phase transition 
at a critical force F c with the oscillator rest state as its 
order parameter. We find that the DNA force-extension 
curves correspond to cycling at different rates the curves 
of the first order phase transition. As in the experiments, 
we find a region of metastability in a certain range of 
forces, close enough to F c . Moreover, the residence time 
spent at the basin of attraction of both the stable and 
metastable states obey an Arrhcnius law: its logarithm 
is proportional to F — F c . 

The rest of the paper is as follows. In Section [TTJ the 
oscillator-spin model is motivated in a biological context 
and its equilibrium properties analyzed. The dynami- 
cal behavior of the model is analyzed in Section |III| The 
oscillator obeys Newton's second law with a mean-field 
force due to the coupling with the spins. The latter flip 
stochastically following Glauber dynamics at tempera- 
ture T . This causes the oscillator position to become 
a stochastic process. In the limit of fast spins, the spin- 
oscillator coupling gives rise to a nonlinear friction term, 



which drives the oscillator to equilibrium. In Section IV 



we present Monte Carlo simulations of the system and 
analyze them in the light of the effective potential acting 
on the oscillator. Section [V] contains the main conclu- 
sions of the present work. 



II. THE MODEL 

We consider a one dimensional chain of length Lq. 
There is a large number N + 1 of "internal" degrees of 
freedom sitting at regularly spaced lattice sites that are 
modeled by Ising spins. Thus, the distance between spins 
is do — Lq/N (see fig. [I]). Assume that we stretch the 
chain so that its length becomes L = Lo+A. For the sake 
of simplicity, we will assume that the spins are regularly 
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FIG. 1. Sketch of the model described in the main text. 
The spin representing the internal degree of freedom at each 
lattice site is shown. 



spaced after the stretching, so that the distance between 
two neighboring spins changes to d — do + A/N. This 
assumption amounts to a "mean field" approximation. 
The potential energy of the system is 

l N+l 

V( A, tr) = -muj 2 A 2 + J( A) a t a l+ll (la) 



where 



J(A) = J - fiA, 



(lb) 



is a function of A. The potential V contains a harmonic 
macroscopic elastic term mujA 2 /2 and a spin energy aris- 
ing from a nearest-neighbor interaction. The spin cou- 
pling constant J depends linearly on the separation be- 
tween sites, and equals Jo f° r the initial chain length L$. 
This simple choice is reasonable for A <C L. The inter- 
action between nearest neighbor spins mimics (in a very 
simple way) the short-ranged interaction between the in- 
ternal degrees of freedom of complex biological molecules 
like nucleic acids. We assume that both J and [i are pos- 
itive. Let us define 



= A- 



Jo 
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(2) 



such that J(x) vanishes for x — 0. For x < (folded 
state) the interaction between the spins is ferromagnetic, 
while for x > (unfolded state) it is antiferromagnetic. 
In terms of a;, V becomes 
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V = -muj 2 x 2 + F c x - [ix >| <Ji(T l+1 , F c = 

1 i=i P 

except for an irrelevant additive constant. The parameter 
F c has the dimensions of a force. If an external load F is 
applied to the system, a new term —Fx is added to ^ 
so that the potential energy is now 

1 



N+l 



muj Jq 



N+l 

V = -mui 2 x 2 — Hx — [ix (jjcr,-_|_i, H = F — F c . (4) 

i=l 



This potential energy is the same as introduced in ref. 
[TUIU2) . except for the extra term —Hx. Interestingly, in 
the "zero field" case, our system has a second order phase 
transition at a critical temperature T c , given by [I0HT2] 



T c 



M 2 (iV+l) 
muj 2 kB 



(5) 



where ks is Boltzmann's constant. The rest position 
x is the order parameter of the transition: for T > T c , 
x = 0, whereas for T < T c there are two equally probable 
equilibrium states with rest positions x — ±xq (xq > 0). 

The potential energy Q has the key ingredients to 
model DNA/RNA behavior in unfolding/refolding exper- 
iments. For T < T c and applied force F < F c , H < 
stabilizes the solution with rest state — xq (folded state), 
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whereas fori* 1 > F c , H > and the stable solution has 
rest state +xq (unfolded state). The main effect of the 
"external field" H is that the system undergoes a first 
order phase transition at H = (F = F c ) 27 . Thus, F c 
is the a critical value of the force F: at any temperature 
T < T c the oscillator rest position as a function of the 
force changes abruptly at F = F c . Furthermore there is 
a region of metastability around F c , as discussed below 
in Section iHAl 

To analyze our model, it is convenient to render its 
equations of motion dimensionless first. Let uj^ 1 be the 
time unit. The elastic and the spin term in the potential 
are of the same order if the scale of the oscillator position 
is [x] = /j,(N + l)/(mw 2 ). The force and the spin term 
are of the same order provided the scale of F is [F] — 
fi(N + 1), and therefore the order of magnitude of the 
potential is [F][x] = fi 2 (N + l) 2 /(mw 2 ) = (N + l)k B T c . 
This scaling is reasonable because makes V extensive. 
For the critical temperature T c to be size-independent, 
we must assume that fj, scales as (N + l) -1 / 2 in the limit 
of large system size [T01 - H2"] . Thus we can define nondi- 
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J H 


















-Hb 



»(N+1) 1 

m.!,)'^ til 



F V 
fx{N + l) (N+l)k B T c 



TABLE I. Nondimensional units and parameters. 

mensional variables according to x* = x/[x], t* = t/[t], 
V* = V/[V],..., where the units [x], [t], [V],... are as 
defined in Table [I] The dimensionless potential is 



V* = 



V 



,.*2 



(N + l)k B T c 



H* 



with 



F* = 



F c 



1) 



-H*x* - 

F* - F*, 

muj 2 Jo 
= n 2 (N + l) 



AT- 



* N+l 
i=l 

(6a) 
(6b) 

kn-l-r 



We will drop the asterisks in the following (so as not to 
clutter our formulas), and from now on every expression 
will be written in terms of the dimensionless variables 
and parameters. 



FIG. 2. (Color online) Plot of H vs z cq for 9 = 0.5. The 
values ±Ht between which there are three equilibrium val- 
ues of the oscillator position are shown. For a given value 
of H, —Hb < H < +Ht, the two locally stable equilibrium 
points xl and xr (green) and the unstable one xu (red) are 
indicated. The two symmetric stable equilibrium points Lxcq* 
corresponding to the zero field case are also shown. The qual- 
itative shape of the curve is the same for all the subcritical 
temperatures 6 < 1. 



sum over the spin variables to obtain the marginal dis- 
tribution probability 

V cq (x) = J2 ^(x, a) = i cxp [—(N + l)V cS (x)/9] , 

(9) 

where Z = 2 N Z. In Eq. (I9j), V c ff is an effective potential 
for the x variable, 



V ff(aj) = y-Ha:-flliicosh(|), (10) 

whose minima will be the stable equilibrium values of x. 
Therefore, 



A. Equilibrium state. Effective potential 



x cq = H + tanh( 2 ^p) , 



(11) 



In equilibrium, the joint probability distribution for 
the oscillator position x and the spin configuration cr = 
{(Ji, . . . , ctat+i} is the canonical distribution which, in 
nondimensional variables, is 

V cq (x, *) = \ cxp [-(N + l)V(x, a)/9] . (8) 

Here Z is a normalization constant. Let us study the 
equilibrium values of the oscillator position x. Then, we 



gives the oscillator rest position a; oq in equilibrium as a 
function of the dimensionless external field H = F — F c 
and temperature 9. 

For H = 0, we recover the model analyzed in refs. 
[TUIU2) . in which a; oq = is always a solution for any 
9. For 9 > 1, it is the only solution, it corresponds to 
a maximum of V eq and is therefore stable. At 9 = 1 
two new stable equilibria corresponding to two different 
maxima of V cq bifurcate from that having a; cq = 0. For 
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(0) 



9 < 1, the positions of these maxima are ±x eq 
9 — > 1~ we have 



As 



(12) 



On the other hand, for — > + , we have x^ — > 1, which 
is the maximum value of x^ > 0. For H = 0, the two 
(positive or negative) equilibrium rest positions ia;^ are 
equiprobable because they correspond to equally deep 



minima of the effective potential (10), which is an even 
function of x for H = 0. 

For H ^ and temperatures below critical, 9 < 1, the 
field term — Hx in eq. ( 10 1 breaks the symmetry between 



x > and x < 0, the effective potential is no longer an 
even function of x. If H < 0, the field term favors the 
negative branch x < 0, since it gives a negative contri- 
bution to the effective potential. Therefore, we expect 
to find the system in the "folded" state x < for low 
values of the loading force, F <C F c . On the contrary, for 
H > 0, the energy of the positive branch x > will be 
lowered by the field term, stabilizing it. Thus, the system 
will be in the "unfolded" state x > for high values of 
the loading force, F ^> F c . 

Differentiating Eq. (11) with respect to H, we find 



dux, 



1 



1 



sech 



V 9 I 



(13) 



The right hand side of this expression is the reciprocal 
of the second derivative of the effective potential, and 
therefore we can rewrite ( 13 ) as 



d H x cq = 1. 



(14) 



Then the effective potential has a local minimum and 
the corresponding equilibrium rest position is stable for 
dHX e q > 0, while the equilibrium rest state is unstable 
for dnx cq < 0. Since dnx cq = 9/(8 — 1) for x cq = 0, the 
zero rest position of the oscillator is stable at 9 > 1 and 
unstable at subcritical temperatures 9 < 1. 

The equilibrium position x cq is a monotonically in- 
creasing function of the applied load F (or the applied 
field H = F — F c ) for supercritical temperatures, 9 > 1. 
There is a unique value of x eq for each value of the ap- 
plied force F if 6 > 1, and it is stable. As shown in Fig. 
[2] and given by ( p~3[ ) , x cq is not monotonic for subcritical 
temperatures 9 < 1: It has a local maximum at —Xb and 
a local minimum at +Xb, with 



cosh* (5) = i. (15) 
For each value of H between — Hb and +i?b, with 



H b = x b - tanh (Jj-^J 



(16) 



there are three possible values of x eq : xjj between — x b 
and Xb is therefore unstable, while the other two, xl < 
and xr > are locally stable. The absolute minimum 



of the potential corresponds to the value of x having the 
same sign as the applied field H. The other local min- 
imum, with sgn(a;) 7^ sgn(iJ), is a metastable state in 
the thermodynamic sense. Then, we expect to find bi- 
stability in the system for \H\ < Hb, i.e. for a given range 
of loadings \F — F c \ < Hb around the "critical" force F c . 



III. DYNAMICS 

The Hamilton equations of motion corresponding to 
the nondimcnsional Hamiltonian function 



P 

H(x,p,(T) = — + V(x,a) 



with potential energy given by Eq. ( |6a[ ), are 
x=p, p = -d x V(x,a), 

so that 

1 N+l 

x = -x + H + N Wi+i- 



(17) 



(18) 



(19) 



According to the stochastic Glauber dynamics, the spins 
flip at a rate [5S] 



Wi(x,tr) = 



1 — <7j (<7i + a i+ i) 



j(x) = tanh 



2.i; 



(20a) 



(20b) 



Here, Wi(x, <r) is the transition rate from configuration 
<x to Ri<r, the same as er except for the sign of the i-th 
spin. Since the oscillator evolution equation ( 19 ) includes 



a term that depends on the stochastically changing spin 
configuration, the oscillator position becomes an stochas- 
tic process. 

The average values of the spin correlations, 



' iV i-\-n 1 



(21) 



satisfy the system of equations 



At 



1 



,n+l 



+Cj-l,n+l + Ci+i^ n _i)) 



(22) 



for n > 1, with the boundary condition C^o = 1- If 
the oscillator position x were time-independent, the spins 
would reach the equilibrium distribution corresponding 
to the constant value x. Then the average spin corre- 
lations (Ci_ n ) = {t&nh(x/9)] n would be independent of 
i due to the spatial translation invariance. Something 
similar occurs in the limit of large system size, N ^> 1. 
Both x and the spin correlations Cj )7J become macro- 
scopic self- averaging variables, i.e. they tend to their 
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respective "macroscopic values" , x and C n (independent 
of i), which coincide with their averages and are the most 
probable values of the corresponding stochastic variables 
[25] , Splitting both x and the correlations Ci >n in their 
corresponding macroscopic and fluctuating parts, 



x + Sx, Ci 



C n + 6Ci. 



(23) 



where Sx and 8C i>n are ©(A^ 1 / 2 ), Eqs. (|19l), \22\ and 



(23) yield the evolution equations, 



H + Ci, 



(24a) 



dC^ 
dt 



-2C n +j{x) (c n _i +C„+i) 



n > 1, 
(24b) 

to be solved with Cq = 1 and appropriate initial condi- 
tions. In Eq. (24b) we have neglected terms or order 
1/N such as ((Sx) 7 ), (Sx5Ci, n ), etc. 

As the spins represent the internal degrees of freedom 
of the molecule, we assume that they evolve rapidly com- 
pared with the time scale of the macroscopic degree of 
freedom modeled by the oscillator position. Thus, the 
dimcnsionless characteristic attempt rate satisfies a>l 
(recall that the unit of time has been chosen as cu^ 1 ). 
Thus, we can solve approximately the system of equa- 
using a power series in the small parameter 
o ' si) a], a procedure akin to the Hilbert method in 



tions ( |24b 

1 eta 

kinetic theory. If x were time-independent, the spin cor- 
relations C n would reach the equilibrium values corre- 
sponding to x, 



C n . 



eq 



eq 



eq 



X 

tanh ( — 



(25) 



in the long time limit. The leading order correction of 
this result is [TU] 



C\ — Ci :C 



dCi, 



eq 



> Cq dt ' 

where r is the spins average relaxation time |30[ I3T 
1 1 + Cil a 



2a (l-^ q ) 2 ' 



(26) 



(27) 



Equation ( 26 ) does not depend on the initial condition 



for C\ which is forgotten after a time much shorter than 
the oscillator natural period. Inserting Eqs. (25)-(27) 
into (24a), we obtain 



d 2 x 1 l + tanh 2 (f) da? 
dt 2 " + 2^9 \ -tanh 2 (f) ~d~t 



H - tanh - = 0, 



(28) 

which can be rewritten in terms of the nondimcnsional 



effective potential ( 10 ) and the friction coefficient 



R(x) 



1 +tanh 2 (f) 
1 -tanh 2 (|)' 



(29) 



d 2 x 
di 2 " 



1 dx 
— — -R ix) — . 

2a6 v ' dt 



(30) 



This approximate evolution equation gives the dynamics 
of the macroscopic value x of the oscillator position for 
fast spins: the nonlinear friction term drives the system 
towards equilibrium, which corresponds to the minima 
of V e ff. Both the "renormalization" of the potential to 
V e ff and the nonlinear friction term are a consequence 
of the coupling between the oscillator and the internal 
degrees of freedom. Equation (30) ceases to hold as 



9 — > + because the spin relaxation time given by ( 27 ) di- 
verges. A detailed discussion on this point can be found 
in refs. [10\ 111). In the remainder of the paper, we will 
restrict ourselves to a temperature range for which the 
spins change rapidly compared to the oscillator motion 
and eq. (30) holds. 



A. Metastability region 

For H = and subcritical temperatures, 9 < 1, the 
effective potential has two equally deep minima at the 
symmetric positions ±x^ of Section |nj and a maximum 
at x = 0. For \H\ < H b (see figure [2] for a qualita- 
tive picture), the effective potential has two minima at 
xn > and xl < and a metastability region appears. 
The globally stable position satisfies XiH > (i = R,L), 
while the other one is a metastable state in the thermo- 
dynamic sense. It must be stressed that this bistability 
is present for all the subcritical temperatures 9 < 1, it is 
not limited to a region near the critical temperature. 

Let us analyze in more detail the situation for weak 
fields. Expansion of Eq. ( 11 ) in powers of H = F — F c <C 
1 gives 



XRjL =±x<g + X H + 0(H 2 ), (31) 
he 

field, and x is the zero-field "susceptibility 



where a^q* is given by the solution of Eq. ([11]) for zero 



X = 9 H x cq \ H=Q = 1 



1 x 



(oy 

eq 



(32) 



20^ to lowest order in H. Close 



Therefore, xr — xl 
to the critical temperature, 9 — > 1~, substitution of Eq. 
([12]) into p2l) yields x ~ 9/[2{l-9)]. On the other hand, 



the unstable position xu changes from zero to 



x u 



1 



1)H + 0(H 2 



(33) 



According to (J9j) , the transitions between the two minima 
are hindered by the presence of large energy barriers, 

Br,l = (N + 1) [Vc ff (^) - V eS (x R , L )} , (34) 
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which are proportional to the system size N + 1. For 
H = 0, both equilibrium states have the same energy 
barrier, 



B^ = -(N+l)V ea (x^)>0. 



(35) 



For H ^ 0, the barrier corresponding to the stable equi- 
librium point verifying XiH > (i — R, L) is larger than 
the one for the metastable equilibrium point. The bar- 
rier from the metastable state tends to zero as \H\ — > Hi, 
because xu and the metastable equilibrium position [xl 
for H > 0, xr for H < 0) coalesce in that limit. The 
maximum effective potential at xjj separates the basins 
of attraction of the minima at and %r- The system 
spends long periods of time oscillating in the vicinity of 
either xl or xr until it is able to hop to the other min- 
imum via a thermally activated process. The residence 
times in each basin of attraction increase exponentially 
with N + 1, and we have to consider a system of moder- 
ate size in order to see hopping between the two minima 
on a reasonable time scale. 

Let us estimate the barrier height for weak fields. From 
Eqs. (10) and (31 1, we have 



"eq 



T Hx 



(o) 

cq 



0(H 2 ), 



(36) 



while V e ff(xu) = 0{H 2 ). In Eq. Q and for the rest 
of this section, the upper sign corresponds to xr and 
the lower sign to Xl- The respective barriers from the 



equilibrium states xr l defined in Eq. (34) are 



B r<l ~bW±(N + 1)HxW. 



(37) 



The residence times in the respective basins of the min- 
ima should have the Arrhenius form, 



tr,l = T exp(B RtL /9), 



(38) 



where tq is some characteristic time. By inserting Eq. 
d37b into p8|, we get 



f, (iv+i)4V 

tr.l = r c exp ± -H 



(39) 



where r c = To exp(B^/9) is the residence time for zero 
field in each basin (the same for them both). Equation 
(39) is the main result of this section; we should stress 
that it is valid for weak fields H -C 1, but (N + l)H 
can be of the order of unity or even a large number. For 
H > 0, we have tr > tl, because xr is the globally 
stable state, while for H < it is tr < tl, since xr is 
the metastable state in that case. Interestingly, the ratio 
of the average lifetimes tr/t l gives the so-called equilib- 
rium constant K for folding/unfolding at that force value 
|16j . Therefore, for our model we arrive at 



rr TR 

K = — = exp 

TL 



'2(N+l)x^ 



H 



H = F-F C1 (40) 



In K is a linear function of the applied force F, A com- 
pletely analogous behavior has been observed experimen- 
tally [IB]. The exponent in Eq. (40) is simply N+l times 



the difference of values of the effective potential (which 
plays the role of the free energy per particle) between 
both states, as readily seen by making use of Eq. (36). 



IV. NUMERICAL RESULTS 

We have performed Monte Carlo simulations of the 
system dynamics introduced above. In all the cases pre- 
sented here, the dimensionless temperature has been cho- 
sen to be 9 = 0.9. The qualitative shape of the equilib- 
rium H vs. x curve is similar to the one shown in Fig. 
|2j For 9 = 0.9, the oscillator rest position at equilibrium 
and zero field is ±xiq ~ ±0.525, as given by eq. ( fl3~| ) 



Interestingly, the approximation in eq. ( 12 1, which is ex- 



pected to be valid very close to the critical temperature, 
gives quite a good estimate x^ ~ 0.547. The points at 
which the curve H vs. x has either a maximum or a min- 
imum are ±Xb = ±0.295; the corresponding values of the 
field are T^b = ±2.15 x 10~ 2 . There is metastability for 
applied fields in the interval \H\ < Hi,. 



A. Force-extension curves 

Let us analyze the force extension curves of the model. 
The pulling/pushing cycle is as follows. We start at 
the equilibrium configuration corresponding to a; m i n < 
(folded state), so that the initial value of the field is 
tanh(a; m i n /0) given by eq. (fTTJ . We pull 
the system by a stepwise increment of the field: at each 
step, the field is increased by AH, and then the system 



is allowed to evolve during a given time At. At the end 
of this period, we record the oscillator position x. Then 
we increase again the field by AH and continue the pro- 
cess in the same vein. The pulling process ends when we 
reach a positive value of the field -ff max = — -ffmin- Then, 
we start to push back, decreasing the field by AH at 
each step, until we reach again the minimum field value 
-ff m i n . As during the pulling process, we record the value 
of x at fixed H after the evolution time At. This process 
is completely analogous to that carried out in unzipping 
experiments with biomolecules. 

Figure [3] shows some typical pulling/pushing cycles 
with different loading rates AH/ At. We have used a 
system with N ± 1 = 1000 spins, and a spin attempt rate 
a = 4, large enough for the spins to be fast as compared 
to the oscillator [TD]. For all the curves, the minimum 
value of the oscillator position is x m { n = —1.5, and the 
field increment at each step is AH = 10~ 3 , which is 
smaller than Ht, and it allows the system to visit the 
metastability region. The system behavior is qualita- 
tively similar for other parameter values, as long as the 
temperature 9 < 1. The loading rate is changed by vary- 
ing the amount of time At at each step. The red (solid) 
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0.15 




FIG. 3. (Color online) Hysteresis cycles for different values 
of the loading rate AH/ At. The red solid lines correspond 
to the unfolding process (At = 5, 100 and 10 4 from top to 
bottom), and the green dashed lines to the refolding process 
(At = 5, 100 and 10 4 from bottom to top). The H versus x 
curve at equilibrium is plotted with blue stars. 



lines correspond to the unfolding process (At = 5, 100 
and 10 4 from top to bottom), and the green (dashed) 
lines to the refolding process (At = 5, 100 and 10 4 from 
bottom to top) . These numerical curves are qualitatively 
similar to those observed in unzipping experiments with 
nucleic acids [TB]. There is always some hysteresis as the 
unfolding and folding curves are not superimposed. The 
area of the hysteresis cycle increases with the loading 
rate, being large for the largest loading rate considered 
and almost zero for the smallest one. The main differ- 
ence with the experimental results is that, in our model, 
the extension of the molecule is not linked to a drop- 
ping of the loading force (in order to see this effect in a 
real experiment, see for instance figs. 2(A) and 2(E) of 
ref. [16]). In the experiments of Ref. [16], the total length 
between the beads localizing the molecule is controlled. 
This corresponds to fixing the length L in our model, not 
the load given by H as we have done. When the force 
is externally controlled in experiments, there is no drop 
of the loading force at the extension transition and the 
hysteresis cycle described by the molecule extension is 
completely analogous to ours [25] . 

Our numerical results can be explained as follows. As 
depicted in Fig. [2] and observed in the last paragraph of 



Sec. II A for \H\ > H b we have a unique stable equilib- 



rium point for \H\ > Hb, the zipped state x cq < —Xb < 
for H < —Hb < and the unzipped state x cq > +Xb > 
for H > Hi,. For \H\ < Hb, there are two (locally) stable 
rest positions xl < and Xr > and an unstable po- 
sition xjj between them. The true thermodynamic equi- 
librium state of the system corresponding to the global 
minimum of the potential satisfies x^H > (i = L, R) 
whereas the metastable state has rest position such that 
XiH < 0. 

Let us consider again the pulling processes in fig. [3] 



starting from the equilibrium configuration correspond- 
ing to a low value of the applied force, H < 0, \H\ 3> \Hb\. 
When we increase the force at a moderate rate, the oscil- 
lator follows the equilibrium curve for negative values of 
H, with x eq < 0, because its relaxation time is small com- 
pared to At, and the friction term in cq. ( |30[ ) can drive it 
to equilibrium. When the field reaches H = (or a very 
small value), the stable equilibrium position of the os- 
cillator changes discontinuously from —xiq to aieq . The 
important question is now whether At issufficiently long 
for the oscillator position to overc ome the energy bar- 
rier B^ at zero field, given by Eq. (35 1. If the answer is 



positive, the system jumps during the time interval At to 
the other stable branch where x > 0, and it stays on that 
branch when the force is further increased. In the push- 
ing back experiment, the system reverses its path and the 
behavior is almost reversible. This behavior is observed 
for the largest value At = 10 4 , corresponding to a loading 
rate AH/ At = 10~ 7 . For larger pulling rates, such as the 
other ones considered in the same figure, the system does 
not have enough time to jump over the energy barrier at 
H = 0. Therefore, x moves over the metastable branch 
with XiH < 0, until the barrier decreases sufficiently for 
the oscillator position to jump to the most stable branch, 
with x > 0. Of course, the actual part of the metastable 
region visited by the system depends on the loading rate. 
For the highest loading rate considered, corresponding to 
At = 5, the system visits the whole metastable branch 
up to the maximum. A similar line of reasoning explains 
the behavior observed in the refolding curve. It is in- 
teresting to note that the hysteresis cycle found for high 
loading rates is not a non-equilibrium behavior, as previ- 
ously suggested [T6HI9] , but it arises from the sampling of 
the regions of metastability for subcritical temperatures. 

A pulling experiment corresponding to a rate even 
slower than the smallest loading rate in fig. [3] is plotted 
in figure [4] Again, the minimum value of the oscillator 
position has been chosen to be a; m j n = —1.5, and the field 
increment at each step AH = 10~ 3 , but the time spent by 
the system at each value of the force is very large, namely 
At = 10 5 . For the sake of clarity, the region around the 
critical force F — F c (H — 0) has been zoomed in. Wc 
observe several jumps between the folded (x < 0) and 
unfolded (x > 0) states for \H\ = \F - F c \ < 5 x 10~ 3 . 
This behavior is a clear signature of the bistability shown 
by our system in the region \H\ < Hb- At zero field, the 
energy barriers between the two stable states are identi- 
cal, Tfj = tl = t c in ( 39 ) , and the reverse transition is 
equally likely. For At 3> r c , the oscillator position has 
enough time to surpass the energy barriers several times 
and it may go back and forth from one state to the other, 
as shown in Fig. [4] 



B. Constant force experiments 

We have also carried out Monte Carlo simulations at 
constant force. We have chosen 8 — 0.9, a small value 
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FIG. 4. (Color online) Detail of the metastability region 
| if | < Ht for a pulling experiment with a very slow loading 
rate, AH = 10 -3 and At = 10 5 . Hopping between the zipped 
and unzipped state is clearly seen for \H\ = \F — F c \ < 5 x 
10~ 3 . 




5000 10000 15000 20000 25000 
t 

FIG. 5. (Color online) Time trace of the oscillator position 
for a constant force experiment in the metastability region, 
namely H = 3 x 1CP 3 . Hopping between the two locally 
stable equilibrium points of the oscillator is observed. 



of the field, H = 3 x 10~ 3 < Hf,, and a smaller size, 
N + 1 = 500, in order to keep the simulation time un- 
der control. Equation ( [TT| ) gives the two locally stable 
oscillator rest positions xl = —0.509 (metastable) and 
Xr = 0.540 (globally stable), separated by the unstable 
oscillator position xjj — —0.027. These values agree with 
the weak field expressions (pi)) -(33). Figure [5] shows a 
time trace of the system. The oscillator jumps stochas- 



tically between the two values xl and xr correspond- 
ing to the two locally stable equilibrium points. The 
system spends more time in the stable state xr (recall 
H = 3 x 10~ 3 > 0), because it has to surpass a larger 
energy barrier in order to escape therefrom. 

In order to check Eqs. (39) or (40) for the residence 
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FIG. 6. (Color online) Logarithm of the ratio K = tr/tl 
(crosses) of the residence times for the two equilibrium states 
xr,l in the metastability region. The number of spins is iV + 
1 = 500. Also plotted is the best fit to the Arrhenius law 



(solid line), Eq. (40) 



times in each basin of attraction, we have measured the 



average time spent in each basin for different values of the 
applied field in the metastability region \H\ < Hi,. Wc 
find that Ihtr^l increases linearly with H, with a slope 
that agrees with the theoretical prediction of Eq. ( 39 1 . 
We have plotted the ratio of average lifetimes K — tr/t^ 
defined in Eq. (40) as a function of the applied field H 
in Fig. H The rein, In K shows a linear behavior, simi- 
larly to that seen in actual experiments [XT] . The slope 
m = dlnAT/d-ff obtained numerically, m = 542.5 agrees 
well with the theoretical prediction calculated from Eq. 
(40), m = 583.8. This result strongly supports the phys- 
ical picture developed in Section |III A| for the behavior 
of the system in the metastability region, both from a 
qualitative and a quantitative point of view. 



V. CONCLUSIONS 

We have modeled DNA folding/unfolding under an ex- 
ternal load force by a macroscopic linear oscillator cou- 
pled with internal degrees of freedom represented by 
Ising spins that undergo Glauber dynamics. The simple 
mean-field character of the model prevents us from doing 
quantitative comparisons with the real experiments, and 
we have to settle for qualitative comparisons. We can- 
not simulate position-controlled experiments, only force- 
controlled ones. 

Despite its limitations, the picture arising from this 
simple model is physically appealing. The hysteresis cy- 
cles show that the system exhibits a metastable equilib- 
rium behavior in the unzipping experiments, not a true 
non-equilibrium behavior, as it was suggested previously 
[T(3lU9| . In this regard, the unfolding/refolding cycles are 
quite different from the truly nonequilibrium hysteresis 
cycles exhibited by glass formers in cooling/heating pro- 
cesses (see [31] and references therein). In the cooling 
process, the glass formers depart from the equilibrium 
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curve and end in a far from equilibrium state at low 
temperatures. In the reheating process, they return to 
equilibrium following a curve different from the cooling 
one, which typically overshoots the equilibrium curve. In 
contrast with this behavior, the hysteresis cycles in our 
system arise because the metastable equilibrium branches 
of the x vs H curve are swept at high loading rate (such 
that the system does not have enough time to find the 
true minimum of the potential in the bistable field inter- 
val | if | < Hb). Due to the shape of the equilibrium x vs 
H curve in fig. |2j the system unzips at a higher value 
of the field than the one at which it rezips. We expect 
that this general picture remains valid for more realistic 
models and/or actual biomolecules. 

As in the experiments, the system hops between the 
zipped and the unzipped state at very small loading rates. 
It has then enough time to surpass the energy barrier 
separating the unfolded and the folded states. When the 
force is held constant with \H\ < (bistable region), 
the system hops back and forth between the two possible 
equilibrium values of the oscillator position. The aver- 
age lifetimes show an Arrhenius-like dependence on the 



applied field H = F — F c , again in agreement with the 
behavior of real systems. 

For all subcritical temperatures T < T c , these behav- 
iors occur due to the first order transition and its asso- 
ciated region of metastability. The phase transition is a 
consequence of the coupling between the oscillator and 
the Id spin system, which introduces an effectively long- 
range interaction between the spins. Similar hidden Id 
long range effective correlations enable phase transitions 
in biological systems [55] , 
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